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The behavior of a family of dissipative dynamical systems representing transformations of two- 
dimensional torus is studied on a discrete lattice and compared with that of conservative hyperbolic 
automorphisms of the torus. Applying dissipative dynamical systems to generation of pseudorandom 
numbers is shown to be advantageous and equidistribution of probabilities for the sequences of bits 
can be achieved. A new algorithm for generating uniform pseudorandom numbers is proposed. The 
theory of the generator, which includes proofs of periodic properties and of statistical independence 
of bits at distances up to logarithm of mesh size, is presented. Extensive statistical testing using 
available test packages demonstrates excellent results, while the speed of the generator is comparable 
to other modern generators. 



I. INTRODUCTION 

Pseudorandom number generation is an important 
component of any stochastic simulations such as molecu- 
lar dynamics and Monte Carlo simulations [l[ . The prob- 
lem of design of reliable and fast generators is of great 
importance and attracts much attention Q. 

There are numerous papers where chaos is considered 
as a requirement for good pseudorandomness. Many 
properties of chaotic dynamical systems are discussed 
in this respect: ergodicity, sensitivity to initial condi- 
tions, mixing property, local divergence of trajectories, 
deterministic dynamics and structural complexity. These 
properties resemble certain properties of pseudorandom- 
ness and are considered in the literature as desirable 
properties for pseudorandomness. Several pseudoran- 
dom number generators based on chaotic maps have been 
proposed in the literature [H, H|. However, the behav- 
ior of dynamical systems on a discrete lattice is stud- 
ied much less than in continuous space and a number of 
corresponding important questions still remain open. In 
this work I show that applying dissipative dynamical sys- 
tems to pseudorandom number generation can result in 
substantially preferable statistical behavior of the corre- 
sponding pseudorandom number sequences, compared to 
applying conservative dynamical systems. 

The present approach extends the method of pseudo- 
random number generation of Ref. [J, [5( , which is based 
on evolution of the ensemble of dynamical systems. Sev- 
eral generalizations are carried out. The connection be- 
tween the statistical properties of a generator and geo- 
metric properties of the corresponding map is uncovered. 
New pseudorandom number generator is proposed. Us- 
ing SSE2 technology, which is supported by all Intel and 
AMD processors fabricated later than in 2003 [f| 0] , ef- 
fective implementations are developed. 

One of the most important properties characterizing 
the quality of pseudorandom sequences of numbers is 
the high-dimensional uniformity and the corresponding 
equidistribution property @ . Unlike other essential char- 



acteristics of pseudorandom number generators such as 
the period length, which is studied in detail in relation- 
ship to nearly all known generators, there are not so many 
examples in which the high-dimensional equidistribution 
property was proved . 

In this paper the proper choice of parameters is estab- 
lished, which results in the validity of the equidistribu- 
tion property for the proposed generator. In particular, 
it is shown that the determinant of the transformation 
has to be an even integer in order for the property to 
hold. The equidistribution is established on length up to 
a characteristic length t: for n < £, each combination of 
successive n bits taken from the RNG output occurs ex- 
actly the same number of times and has a corresponding 
probability 1/2™. The length £ turns out to depend lin- 
early on t, where the mesh size g (i.e. the modulus of the 
basic recurrence) is equal to p ■ 2 and p is an odd prime. 
In other words, for given p, one has £ cx log g. Numerical 
results show that the equidistribution property still ap- 
proximately holds with high accuracy beyond the region 
of its strict validity under the condition n < 6.81ogp. 

I have constructed several realizations for the proposed 
generator (see Table [I]). It is shown in Proposition 2 in 
the section on geometric and statistical properties that 
for the realizations either £ = 2t — 1 or ^ = (t — l)/2 
takes place. The speed and statistical properties of the 
constructed generators are compared with those of other 
modern generators (see Tables HH IIIII) . Practically, the 
generators with smaller values of t (e.g. with prime g) 
also have very good properties for a particular choice of 
parameters, while the generator period is not less than 
p 2 — 1 and increases significantly with increasing p. For 
this reason two realizations with small t are also thor- 
oughly tested. 

Among several statistical test suites available in the 
literature, TestUOl is known to contain very stringent 
batteries of tests for empirical testing of pseudoran- 
dom numbers. At present there are only several known 
pseudorandom number generators that pass all the tests 
even in the sense that no p- value is outside the interval 
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[1CT 10 , 1 - 1CT 10 ] [Hj]. Statistical testing with TestUOl 
confirms excellent statistical properties of the proposed 
realizations. 

The results obtained have further perspectives in view 
of generating large number of guaranteed statistically in- 
dependent pseudorandom streams, which can be partic- 
ularly well-suited for use in a parallel, distributed envi- 
ronment. 



II. THE GENERATOR, ITS INITIALIZATION 
AND PERIOD 

It is suggested in [H, 0] to construct RNGs based on 
an ensemble of sequences generated by multiple recur- 
sive method. The state of the generator consists of the 
values xi^jXj 2 ' <E {0, 1, . . . ,g— 1}, % = 0, 1, . . . , s — 1. 
The transition function of the generator is defined by the 
recurrence relation 

x[ n ^ = kxf 1 ^ — qxf 1 2 "* (mod g), (1) 

where i — 0, l,...,s — 1. The values \ i = 
0,1, ... ,s — 1 can be considered as ^-coordinates of s 
points (x\ n \ y!- ) T , i = 0, 1, . . . , s — 1 of the g x g lattice 
on the two-dimensional torus, then each recurrence rela- 
tion describes the dynamics of x-coordinate of a point on 
the two-dimensional torus: 

fx {n) \ fx {n ~ 1] \ 

Lh)= M [ ? >-i)J (mod3) ' (2) 

where matrix M — ( mi m2 ) is a matrix with in- 
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teger elements, k = Tr M, q = det M and Tr M 

is a trace of matrix M [1, [l3|, [l4j. Indeed, it 

follows from ([2]) that kx\ n ^^ — qx[ n ~ 2 ^ = (mi + 

\ (n-l) , \ (n-2) , (n) (n-l)x , 

m^x} -(mim 4 — m 2 m 3 )x i = \%l -m 2 yl ) + 

(n-1) (n-2) (n-2) (n) 

I (n-l) (n-2).. , (n-1) (n-2) , (n) 

m^m^yf 1 ^ + ra^m^yi 1 ^ = x\ n ^ (mod g). The ba- 
sic recurrence ([T]) is therefore closely related to so- 
called matrix _generator of pseudorandom numbers stud- 
ied in @, [H, [l| . 

The output function is defined as follows: 

a (»>=2N n) /sJ-2\ (3) 

i=0 

where i — 0, l,...,s — 1, i.e. each bit of the output 
corresponds to its own recurrence, and s = 32 recurrences 
are calculated in parallel. 

For g = p ■ 2 , where p is a prime number, the char- 
acteristic polynomial f(x) = x 2 — kx + q is chosen to be 
primitive over Z p . Primitivity of the characteristic poly- 
nomial guarantees maximal possible period p 2 — 1 of the 
output sequence for g = p. It is straightforward to prove 
that taking g = p ■ 2* instead of g = p does not reduce 
the value of the period. 



TABLE I: Parameters of the new generators. 



Generator 


9 


k 


<1 


V 


Period 


GM29.1-SSE 
GM55.4-SSE 
GQ58.1-SSE 
GQ58.3-SSE 
GQ58.4-SSE 


2 29 -3 
16(2 51 - 129) 
2 29 (2 29 - 3) 
2 29 (2 29 - 3) 
2 29 (2 29 - 3) 


4 

25( 
8 
8 
8 


2 

176 

48 
48 
48 


1 

4 
1 
3 
4 


= 2.8 • 10 17 

> 5.1 ■ 10 30 

> 2.8 ■ 10 17 

> 2.8 • 10 17 

> 2.8 ■ 10 17 



There is an easy algorithm to calculate x^ in (fTJ) 
very quickly from a;^ ) and x^- 1 ' for any large n. In- 
deed, if x^ = k n X^ - q n x^ (mod g), then x^ = 
(k 2 — 2q r Ax^ 2n ^ — q 2 x^ (mod g). As was mentioned al- 
ready in Q, this helps to initialize the generator. To ini- 
tialize all s recurrences, the following initial conditions 

are used: xf ] = x^ A \ x[ 1] = x^ A+1 \ i = 0, 1, . . . , s - 1. 
Here A is a value of the order of (p 2 — l)/s. The author 
has tested realizations with various values of A of the 
order of (p 2 — l)/s and found in all cases that the spe- 
cific choice of A was not of importance for the properties 
studied in the next sections. Short cycles and, in partic- 
ular, the cycle consisting of zeroes, are avoided if at least 
one of x^ and x^ is not divisible by p. As a result of 
the initialization, all s initial points belong to the same 
orbit on the torus of the period p 2 — 1, while the minimal 
distance A between the initial points along the orbit is 
chosen to be very large. 

In addition to the realizations based on the output 
function (J3]) that takes a single bit from each linear re- 
currence, I have also constructed realizations based on a 
more general output function 

a<*> = Y j \_2"xf ) /g\-2™, (4) 

i=0 

where v bits are taken from each recurrence and i — 
0, 1, . . . , s — 1. For example, GM55.4-SSE realization cal- 
culates only s = 8 recurrent relations in parallel and takes 
v = 4 bits from each number. Pseudorandom 32-bit num- 
bers can be generated if sv > 32. The sequence of bits 
{ [2 v x\ n '/g\ }, where i is fixed and {x^} is generated 
with relation @ will be designated below as a stream 
of v-bit blocks generated with matrix M. The pairs 
x[° , Xi G TL g for the recurrence ([1]) and x\ , y 8 - € Z g 
for the recurrence represent seeds for the streams of 
u-bit blocks generated with (fTJ and ([2]) respectively. Con- 
sider the set of admissible seeds containing all seeds such 
that at least one of the two values is not divisible by p. 
Selecting the seed at random from a uniform distribution 
over the set of admissible seeds determines the probabil- 
ity measure for output subsequences of a stream of v-bit 
blocks. Such probabilities are considered below in the 
next section. 

The parameters for the particular constructed realiza- 
tions of the generator are shown in Table H The parame- 
ters are chosen in order for the characteristic polynomial 
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x 2 — kx + q to be primitive over Z p . In addition, as is 
shown below, value of q must be divisible by 2 V in order 
for the equidistribution property to hold. Also the value 
of (k + q)g should not exceed either 2 32 or 2 64 in order to 
effectively calculate four 32-bit recurrences or two 64-bit 
recurrences in parallel within SIMD arithmetic. In the 
particular case t = and v = 1 the method reduces to 
that studied earlier in [J, |5| . Program codes for the new 
generators and proper initializations are available in [l6| . 



III. GEOMETRIC PROPERTIES AND 
STATISTICAL PROPERTIES 

In 0] a connection is established between statistical 
properties, the results of a random walk test and geo- 
metric properties of the cat maps. Cat maps are simple 
chaotic dynamical systems that correspond to transfor- 
mations ^ for q = dctM = 1, i.e. hyperbolic auto- 
morphisms of the two-dimensional torus. In particular, 
it is proved in Q that the probability of sequence 0000 
of the first bits generated by a single cat map depends 
only on the trace k of a matrix M and for even k is 
equal to P = P k 2 /(k 2 - 1), where P = 1/16. If k 
is odd, then all sequences of length 4 are equiprobable. 
The probability of sequence 00000 of length 5 is equal to 
P = P (l + l/(3fc 2 - 6)) for odd k, where P = 1/32. 
The condition P > Pq signifies that the 5-dimensional 
equidistribution never takes place for q = 1, i.e. for con- 
servative hyperbolic automorphisms of the torus. In this 
work a more general case q ^ 1 involving dissipative dy- 
namical systems is studied. 

Fig- HI shows the regions on the torus obtained in |4| for 
the third points of sequences of length 5 for the matrix 
( 1 2 ). The regions correspond to the sequences of length 
5 of the first bits generated by the respective RNG, and 
the areas of the regions are equal to the probabilities of 
the sequences. Each region is drawn with its own color. 

Let X, = {(x,y) T \i/2 v < x/g < (i + l)/2 v ,0 < 
y/g < 1}, i.e. the torus is divided into 2" vertical stripes 
Xq,Xx, . . . ,X2 v -i of equal area. Suppose that g is di- 
visible by 2 V and consider the shift S : (x,y) T — > (x + 
g/2 v ,y) T (mod g), i.e. S(Xi) = A (4+1) ( mod 2 «). The 
shift S is a superposition of two rotations: S = R2R1, 
where i?i is a 180-degree rotation with respect to the 
point (l/2" +1 , 1/2) T and R 2 a 180-degree rotation with 
respect to the point (l/2 v , 1/2) T . 

Proposition 1. If (1) M = ( mi m2 ) is a matrix with 
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integer values mi, m,2, m.3, m^, (ii) mi, q = detM and 
g are divisible by 2 V , (iii) the image of the lattice g x g 
with the transformation M J is invariant with respect to 
the shift S for j = 0,1, ... ,n, then all the sequences of 
length n in a stream of u-bit blocks generated with matrix 
M are equiprobable. 

Proof. In this case the element m|[ of matrix 

((n) (n) x 
„) („) (mod g) (5) 



FIG. 1: (Color online) The regions on the torus obtained 
in for the third points of sequences of length 5 for the 
matrix (I ^J. Coordinates x/g,y/g are used. These regions 
correspond to the sequences of length 5 of the first bits gen- 
erated by the corresponding RNG. Each region is drawn with 
its own color. 




qm\ (mod g). Hence mj[ is divisible by 2 V for any 
integer n > 1. 

Since is divisible by 2", one has 

M n S{x,y) T = M n (x + g/2 v inod g),y) T = 

M n {x,y) T + (0 1 m^ l) g/2 v ) T . Hence, the set of points A 
such that A £ Xi and M n (A) € Xj passes with the shift 
S into the set of points A such that A S Xu+i} t mo d 2 ") 
and M n {A) G Xj. 

Let's now prove by induction that all sequences of 
length n are equiprobable. Obviously, if g is divisible 
by 2 V , sequences of length 1 are equiprobable: -P(O) = 
P(l) = ... = P(2" - 1) = 1/2". Assume that all se- 
quences of length n — 1 are equiprobable. Let cti = 
P{ix\ . . . x n -i), i = 0, 1, . . . , 2" — 1 be probabilities of se- 
quences of length n. Then ccj = a^+i, i = 0, 1, . . . , 2" — 2 
because the set of points A of the lattice g x g such that 
A e Xi, M (A) eX Xl ,..., M n ' x {A) e X Xn _, passes with 
the shift S into the set of points A of the lattice gxg such 
that A g X (l+1) ( mod 2 „), M{A) G X X1 ,. . . , M n ~ 1 (A) G 

-^x n _i- On. the other hand, X^=o 1 ai 1S ^ ne probability of 
sequence x\ . . . x n -\ of length n — 1 and equals 1/2"(™ _1 ). 
Therefore, a, = l/2 vn , i = 0, 1, . . . , 2 V - 1, and all the se- 
quences of length n are equiprobable. Proposition 1 is 
proved. 

The condition that the image of the lattice gxg with 
the transformation M J is invariant with respect to the 
shift S for j = 0, 1, . . . , n, is used in the above consider- 
ation and is necessary for the Proposition 1. For j = 
the invariance means that g is divisible by 2". If g and 
mi are divisible by 2" , then the number of points A of 
the lattice gxg such that A G X and M n (A) £ X is 
equal to the number of points A of the same lattice such 
that A G Xi and M n (A) G X . If g is not divisible by 



4 



FIG. 2: The set of points A such that A £ Xo and M 2 (A) £ Xo (left panel) and the set of points A such that A £ Xi and 
M 2 (A) £ Xo (right panel) for M = ( 22 A and v = 1. Coordinates x/g,y/g are used. 




2" then these numbers are approximately equal because 
the corresponding areas are equal and g is large number, 
and the exact equality holds only if g is divisible by 2 V . 
Fig.Hshows the sets of points {A\A £ X , M 2 (A) £ X Q } 
and {A\A £ X 1 ,M 2 (A) £ X } for M = ( 2 2 2 ) and v = 1. 

Proposition 2. For M = Af = (j 1 " " 2 ) and 



M = 



f244 43\ 
.32 12 , 



the sequences of length 1,2, in a 

stream of bits generated with matrix M are equiprob- 
able, where £ = 2t - 1, £ = (t - l)/2 and £ = (t - l)/2 
respectively. Here 5 = p-2*, where p is an odd prime, and 
the matrices correspond to the realizations GM29-SSE, 
GM58-SSE and GM55-SSE respectively. 

Proof. Let's check that the image of the lattice g x g 
with the transformation M J is invariant with respect to 
the shift for j = 0,1,..., n and n < £. In particular, 
the invariance takes place if there are integers r,l < t 
such that the distance between integer vectors (x + 
g/2 r+1 , y + g/2 l+1 ) T and (x,y) T after applying transfor- 
mation AP is equal to (g/2, 0) T modulo g. This results in 



(m*f> /2 r + mi> /2 l ,m ( 3 3> /2 r + mf /2 l ) T = (1,0) T (mod 
2). For the matrix M = L 2 ) the condition is satisfied 
when r = j/2, I = j/2 — 1 for even j and r = (J — l)/2, 
I = (j + l)/2 for odd j. Thus i = Jrnax + 1 = 2t - 1. 
Similarly, for each of the matrices M = (j 1 ^ n 2 ) and 

M = ( 2 3 4 2 4 1 4 2 3 ) the condition is satisfied for £ = (t - l)/2. 
Proposition 2 is proved. 

Generally, the following statements are also valid. 
Consider a matrix M with integer elements and the fol- 
lowing integer quantities: g — p ■ 2*, q = det M — 
2 u w (mod g), k = TrM = 2 m r (mod g), u > 1, t > v, 
m > 0. Here w, r are odd integers and p is an odd prime. 
Then (i) all 2 3 sequences of length j in a stream of v-bit 
blocks generated with recurrent relation ([T]) are equiprob- 
ablc for j = 1,2,..., L Here £ = \(t - v)/\u/2]] for 
u < 2m and I = \(t — v)/(u — rnj] for u > 2m; (ii) if 
k is even, then the image of the lattice g x g with the 
transformation M 2t is the lattice p x p on the torus; (iii) 
if k is odd, then the image of the lattice g x g with the 
transformation Ml*/"! is not invariant with respect to 
the shift S. 

Although the exact equidistribution property does not 
hold when distance between some points of the sequence 
> 2t 1 numerical results show that the equidistribution 



holds approximately with high accuracy for the sequences 
of bits of length n, where n < 6.8 \ogp. Also, one can take 
n points with arbitrary distances (not exceeding p 2 — 
1) between them along the orbit (i.e. not necessarily 
successive points of the orbit), where n < 6.81ogp, and 
still the approximate equidistribution will hold with a 
high accuracy. The output value in (J3|) consists of 
high-order bits of s = 32 successive points along the orbit 
of matrix M A , where A is the value of the order of (p 2 — 
l)/s. Therefore, according to the numerical results, the 
output value has a uniform distribution with a very 
high accuracy. 

In most cases the image of the lattice g x g on the 
torus with Af- 7 where j > 2t is the p x p-lattice, there- 
fore it is most interesting to study the deviations from 
the equidistribution for the p x p-lattice. I have cal- 
culated the exact areas on the torus which correspond 
to each of the sequences for M — (J 3). The calcula- 
tions were carried out on a PC using Class Library for 
Numbers [I?} for exact rational arithmetics. For each 
of the 2™ sequences of length n = 1,2,..., the corre- 
sponding set of points on the unit two-dimensional torus 
consists of filled polygons. Exact rational coordinates of 
all the vertices of each filled polygon were found. Also, 
the exact number of points of the p x p lattice inside 
each polygon was calculated. The total area of the poly- 
gons for each of the 2 n sequences of length n was found 
to equal l/(2 n ). Such equality of the areas for different 
sequences of the same length was observed for matrices 
with even determinant and was not observed for matri- 
ces with odd determinant. Let A n $, A n> \, . . . , A nt 2^-x 
be the numbers of points of the p x p-lattice correspond- 
ing to the sequences of length n. Then J2i=o 1 = 
p 2 . Therefore, if A n is the set of numbers A n — 
{2"A„, /p 2 ,2"A„, 1 /p 2 ,...,2"A n , 2n . 1 /j ) 2 }, then (A n ) = 
1, where (A n ) is the average value of A n . The dependence 
of logarithm of variance of A n on n is shown in Fig. [3] for 
p = 2 29 — 3. The calculations for smaller values of p and 
larger values of n demonstrate that the dependence of 
log(cr 2 ) on n is almost linear. Calculations show that the 
deviations from equidistribution are negligibly small in 
the sence that a(A n ) is much smaller than (A n ) = 1, for 
n < 6.8 logp. In particular, for p — 2 29 — 3 the deviations 
are small for n < 130. 
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FIG. 3: Variance of the numbers of points of the p x p-lattice 
corresponding to sequences of length n versus n. The values 
are normalized such that (A n ) — 1. 

The variance for the several points of the orbit of ma- 
trix M on the p x p-lattice on the torus, is found to 
substantially depend on the number of points and on the 
value of p, and only weakly depend (within several per- 
cent) on the distances between the points along the orbit. 

IV. STATISTICAL TESTING 

Table ILT1 shows the results of applying the SmallCrush, 
PseudoDiehard, Crush and BigCrush batteries of tests 
taken from [18]], to the generators introduced in Ta- 
ble HI Batteries SmallCrush, PseudoDiehard, Crush and 
BigCrush contain 15, 126, 144 and 160 statistical tests 
respectively. For each battery of tests, Table UU displays 
three characteristics: the number of statistical tests with 
p- values outside the interval [10 -3 , 1 — 10~ 3 ], number of 
tests with p- values outside the interval [10 -5 , 1 — 10 -5 ], 
and number of tests with p- values outside the interval 
[KT 10 , 1 - 10~ 10 ]. Table HI] also contains the results of 
statistical tests for Mersenne Twister generator of Mat- 
sumoto and Nishimira [9j, combined Tausworthe gener- 
ator of L'Ecuyer [ll[ and combined multiple recursive 
generator proposed in [lj|. These generators are mod- 
ern examples of fast RNG implementations with good 
statistical properties (see Sec. 4.5.4 and Sec. 4.6.1 
in (U). Both LFSR113 and MT19937 fail the test 
scomp_LinearComp that is a linear complexity test for the 
binary sequences (see (3), because the bits of LFSR113 
and MT19937 have a linear structure by construction. 
Also LFSR113 fails the test smarsa_MatrixRank (see 
pl|). The period lengths for the generators MRG32K3A, 
LFSR113 and MT19937 are 3.1 • 10 57 , 1.0 • 10 34 and 
4.3 • 10 6001 respectively. 

The usefulness of a RNG for a specific application in 
physics depends on, possibly dangerous interferences of 



the correlations in the specific problem and those of the 
RNG. Modern statistical test suites contain tests that 
reveal known types of correlations for the RNGs, in par- 
ticular, the types that are known to result in systematic 
errors in Monte- Carlo simulations and that were studied 
in 21]. One concludes that the new realizations described 
in this paper possess excellent statistical properties. 

TABLE II: Numbers of failed tests for the batteries of tests 
SmallCrush, Crush, BigCrush [18], and Diehard 0. Testing 
was performed with package TestUOl version TestU01-1.2.3. 
For each battery of tests, three numbers are displayed: the 
number of statistical tests with p-values outside the interval 
[10~ 3 , 1 — 10~ 3 ], number of tests with p-values outside the 
interval [10 -,J , 1 — 1CP J ], and number of tests with p-values 
outside the interval [1CT 10 , 1 - 10" 10 ]. 



Generator 


SmallCrush 


Diehard 


Crush 


BigCrush 


MRG32k3a 


0,0,0 


0,0,0 


0,0,0 


0,0,0 


LFSR113 


0,0,0 


1,0,0 


6,6,6 


6,6,6 


MT19937 


0,0,0 


0,0,0 


2,2,2 


2,2,2 


GM29.1-SSE 


0,0,0 


0,0,0 


0,0,0 


0,0,0 


GM55.4-SSE 


0,0,0 


0,0,0 


0,0,0 


0,0,0 


GQ58.1-SSE 


0,0,0 


0,0,0 


0,0,0 


0,0,0 


GQ58.3-SSE 


0,0,0 


0,0,0 


0,0,0 


0,0,0 


GQ58.4-SSE 


0,0,0 


0,0,0 


0,0,0 


0,0,0 



V. SPEED OF THE GENERATOR 

I have tested the CPU times needed for generating 10 9 
random numbers. The results are shown in Table ILTll for 
Intel Core i7-940 and AMD Turion X2 RM-70 processors 
respectively. The results are presented for different com- 
pilers and optimization options. The compilers in use 
are GNU C compiler gcc version 4.3.3 and Intel C com- 
piler ice version 11.0. The CPU times for the realiza- 
tions GM29.1-SSE, GM55.4-SSE, GQ58.1-SSE, GQ58.3- 
SSE and GQ58.4-SSE introduced in TableUare compared 
with those for Mersenne Twister generator of Matsumoto 
and Nishimira @, combined Tausworthe generator of 
L'Ecuyer [111 and combined multiple recursive generator 
proposed in [19j . 
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TABLE III: CPU time (sec.) for generating 10 9 random numbers. Processors: Intel Core i7-940 and AMD Turion X2 RM-70. 
Compilers: gcc 4.3.3, ice 11.0. 



Intel Core i7-940 


gcc -OO 


gcc -Ol 


gcc -02 


gcc -03 


ice -O0 


ice -Ol 


ice -02 


ice -03 


Source 


MT19937 


13.7 


5.7 


6.9 


2.6 


17.5 


6.5 


2.9 


2.9 


M 


MT19937-SSE 


5.2 


4.8 


5.5 


2.0 


4.9 


4.7 


2.4 


2.0 


\5] 


LFSR113 


10.4 


4.8 


6.8 


3.1 


10.2 


5.0 


4.6 


4.5 


[11] 


LFSR113-SSE 


8.0 


6.8 


6.8 


6.9 


7.3 


6.9 


6.6 


6.5 


[5] 


MRG32k3a 


47.9 


36.3 


35.3 


25.0 


56.1 


33.1 


22.8 


28.1 


[19] 


MRG32k3a-SSE 


9.1 


7.4 


5.8 


5.8 


8.8 


7.4 


6.0 


5.9 


M 


GM29.1-SSE 


22.6 


19.6 


17.5 


18.1 


21.2 


18.7 


18.2 


18.1 


[16] 


GM55.4-SSE 


18.0 


16.8 


15.4 


15.4 


17.7 


16.3 


15.8 


15.7 


[16] 


GQ58.1-SSE 


50.5 


49.2 


47.4 


47.3 


50.5 


48.1 


48.0 


47.7 


[16] 


GQ58.3-SSE 


22.0 


21.2 


19.0 


20.1 


22.5 


20.4 


19.5 


19.5 


[16] 


GQ58.4-SSE 


16.1 


14.7 


12.8 


13.8 


15.5 


13.9 


13.3 


13.3 


[16] 


AMD Turion X2 RM-70 


gcc -00 


gcc -Ol 


gcc -02 


gcc -03 


ice -O0 


ice -Ol 


ice -02 


ice -03 


Source 


MT19937 


31.0 


17.8 


10.8 


7.1 


31.0 


18.7 


5.2 


4.9 


M 


MT19937-SSE 


11.3 


10.3 


11.1 


6.6 


10.8 


9.9 


6.0 


6.0 


[5] 


LFSR113 


14.6 


8.7 


9.6 


5.3 


14.9 


9.1 


6.9 


6.8 




MRG32k3a 


89.0 


60.9 


60.9 


47.0 


89.1 


69.2 


41.5 


41.6 


[19] 


MRG32k3a-SSE 


25.9 


22.3 


18.4 


18.3 


25.6 


22.3 


19.0 


19.0 


[5] 


GM29.1-SSE 


68.5 


64.4 


60.7 


60.7 


67.8 


63.1 


61.7 


61.7 




GM55.4-SSE 


59.8 


54.8 


53.1 


53.0 


58.2 


53.6 


52.8 


52.8 


[16] 


GQ58.1-SSE 


179.6 


179.6 


178.3 


177.8 


183.1 


178.3 


178.5 


178.5 


[16] 


GQ58.3-SSE 


75.5 


73.9 


70.6 


71.1 


74.2 


71.9 


70.4 


70.1 


[16] 


GQ58.4-SSE 


51.9 


51.0 


48.2 


48.1 


53.1 


49.4 


48.2 


48.1 


[16] 



A.R. Bizzarri, J. Phys.: Cond. Mat. 16 (2004) R83. 
[2] D.E. Knuth, The Art of Computer Programming, Vol. 2 

(Addison- Wesley, Reading, Mass., 3rd edition, 1997). 
[3] M. Falcioni et. aL, Phys. Rev. E 72, 016220 (2005); R. 

Lozi, Ind. J. Ind. Appl. Math. 1 (1), 1 (2008); L. Kocarev, 

G. Jakimoski, IEEE Trans. Circ. Syst. 50(1), 123 (2003); 

P. Li et. al., Phys. Lett. A 349, 467 (2006); V. Patidar, 

K.K. Sud, EJTP 6 (20), 327 (2009); N.K. Pareek et. al., 

Int. J. Netw. Sec. 10 (1), 32 (2010). 
[4] L. Barash, L.N. Shchur, Phys. Rev. E 73 (2006), 036701. 
[5] L.Yu. Barash and L.N. Shchur, Comput. Phys. Commun. 

182 (2011), 1518. 
[6] http : //www. intel . com/ support /processor s/pentium4 

/sb/CS-029967.htm 
[7] http : / / support . amd. com/us/Processor_TechDocs 

/24592.pdf 

[8] J.P.R. Tootil, W.D. Robinson, D.J. Eagle, J. ACM 20, 
3, 469 (1973); M. Fushimi, S. Tezuka, Commun. ACM 4, 
516 (1983); R. Couture, P. L'Ecuyer, S. Tezuka, Math. 
Comput. 60, 749 (1993); S. Tezuka, P. L'Ecuyer, ACM 
Trans. Model. Comput. Simul. 1, 99 (1991); P. L'Ecuyer, 
Math. Comput. 65, 203 (1996). 
[9] M. Matsumoto and T. Nishimura, ACM Trans, on Mod. 
and Comp. Sim., 8 (1998) 3. 

[10] F. Panneton, P. L'Ecuyer, M. Matsumoto, ACM Trans. 
Mathem. Software 32(1) (2006) 1. 

[11] P. L'Ecuyer, Math, of Comp., 68 (1999) 261. 



[12] P. L'Ecuyer, R. Simard, ACM Trans. Mathem. Software, 

33(4) (2007) article 22. 
[13] H. Grothe, Statistical Papers, 28 (1987) 233. 
[14] P. L'ecuyer, Comm. of the ACM, 33(10) (1990) 85. 
[15] H. Niederreiter, in Monte Carlo and Quasi-Monte Carlo 

Methods in Scientific Computing, ed. H. Niederreiter 

and P. J.-S. Shiue, Lecture Notes in Statistics, Vol. 106 

(Springer- Verlag, 1995). 
[16] http : //www . comphys . ru/barash/rng_sse2 . zip 
[17] http://www.ginac.de/CLN/ 

[18] P. L'Ecuyer, R. Simard, TestUOl: A Software Li- 
brary in ANSI C for Empirical Testing of Random 
Number Generators (2002), Software user's guide, 
http : //www . iro .umontreal . ca/~simardr/testu01 
/tuOl .html 

[19] P. L'Ecuyer, Oper. Res., 47 (1999) 159. 

[20] P. L'Ecuyer, Chapter 4 of the Handbook of Simulation, 
Jerry Banks Ed., Wiley, 1998, pp. 93-137. 

[21] A.M. Ferrenberg, D.P. Landau, Y.J. Wong, Phys. Rev. 
Lett. 69 (1992) 3382; L.N. Shchur, Comp. Phys. Comm. 
121-122 (1999) 83; P. Grassberger, Phys. Lett. 181 
(1993) 43; L.N. Shchur, J.R. Heringa, H.W.J. Blote, 
Physica A 241 (1997) 579; L.N. Shchur, H.W.J. Blote, 
Phys. Rev. E 55 (1997) R4905; F. Schmid, N.B. Wilding, 
Int.J.Mod.Phys. C 6 (1995) 781. 



